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ABSTRACT 
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I. INTRODUCTION 

Heavy quark effective theories have led to considerable insight into the properties of 
hadrons containing b- and c-quarks. Key among the assumptions underlying the heavy quark 
theory is that one can expand matrix elements in A/ra, where A is some QCD related scale, 
m is the heavy quark mass, and where coefficients in the expansion are independent of m. We 
show in this paper that a similar expansion of the moments of the heavy quark fragmentation 
function leads to valuable information about the heavy quark fragmentation function both in 
leading twist and beyond. We derive a parameterization of the fragmentation function which 
involves heavy quark mass independent parameters and show how the leading parameters can 
be readily extracted from the lowest moments. This permits tests of heavy quark symmetry. 
In particular, the fragmentation function of the b- and c-quarks, both of which have been 
already studied and will be better measured in the future, 1 ' 2 are readily related. We show 
that the commonly used Peterson fragmentation function 3 has a misleading dependence on 
the QCD scale and heavy quark mass; the width is linear and not quadratic in the ratio. 

The basic idea underlying our analysis is very simple. Fragmentation and distribution 
functions measured at a scale Q 2 evolve from "boundary data" — forward matrix elements 
of bilocal heavy quark operators — supplied at some reference scale \i 2 . This boundary data 
incorporates the non-perturbative features of the fragmentation function (at least in regions 
sufficiently far from z = and z = 1 where we will see heavy quark methods and perturba- 
tion theory respectively break down). By evaluating the matrix elements at sufficiently low 
renormalization mass scales, one can exploit heavy quark symmetry to derive a heavy quark 
independent parameterization of the structure function. At this low mass scale, the boundary 
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data can be organized as a power series in A/m, with mass independent coefficients. The mea- 
sured fragmentation function is then obtained by evolving the function to Q 2 via perturbative 
QCD. 

Previous work on heavy quark fragmentation functions either has been completely 
phenomenological 3 or has concentrated on the fragmentation function generated by perturba- 
tive evolution alone. 4 Some authors have explored the possibility of describing the fragmen- 
tation function solely in terms of perturbative QCD. 5 We abandon the objective of previous 
authors of deriving the fragmentation function. We distinguish between perturbative contri- 
butions which can be reliably calculated at high momenta and non-perturbative contributions 
which we parameterize (but don't calculate) by exploiting the heavy quark mass expansion. 
We show that the expansion of the moments of the fragmentation function at any given order 
in A/m only involves a small number of mass independent parameters. Because the same 
parameters enter the b and c quark fragmentation functions and moments of different order, 
tests of heavy quark symmetry will be possible. 

The principle results of our analysis are: 

• At leading twist, we find a simple parametric form for the quark mass dependence of the 
fragmentation and distribution functions at the heavy quark mass scale which applies 
at leading order in A/m. Our form differs from the commonly used function of Ref. [3] 
because our prediction for the width is order A/m, and not order A 2 /m 2 . This difference 
is critical to extracting the scale A which enters the heavy quark expansion. 

• We show that the moments of the fragmentation function can be interpreted in terms of 
mass-independent parameters. We concentrate on only a finite number of these parame- 
ters which we expect can be reliably extracted. This will allow for tests of heavy quark 
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symmetry. Moreover, it could provide further data on which scale A works best in the 
heavy quark expansion. 

• Within our framework it is possible to make rather strong statements about those higher- 
twist effects proportional to powers of m 2 /Q 2 , where Q 2 is the center-of-mass energy 
squared. Since m 5>> A, these are the most important higher-twist corrections. We show 
that the sum of all corrections of order (m 2 /Q 2 ) k to the process e + e~ — > H(P) + X 
vanishes for all k > at leading order in a a . Radiative corrections modify this result and 
give rise to corrections of order -a^Q 2 )^ In We show that these corrections can 
be related back to the leading twist fragmentation function. Mass suppressed corrections 
are of order Am/Q 2 and can also be related back to the leading twist fragmentation 
function. 

The paper is organized as follows. In Section II we derive the leading twist results, first 
for distribution functions, which are pedagogically simpler, although not the case of greatest 
experimental interest, and then for fragmentation functions. In Section III, we investigate 
the applications of our approach. We conclude that a finite number of non-perturbative 
parameters can be reliably extracted by measuring the moments of the fragmentation function, 
using a simple phenomenological form whose parameters are readily interpreted in a heavy 
quark mass expansion. In Section IV we consider higher-twist contributions to e + e~ — > 
H(P) + X, after which we conclude. 
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II. HEAVY QUARK DISTRIBUTION AND FRAGMENTATION 
FUNCTIONS AT LEADING TWIST 

II. 1 Distribution Functions 

Figure 1 shows the virtual Compton scattering diagram for an h quark current whose 
imaginary part yields the dominant contribution to a hypothetical deep inelastic scattering 
experiment from an i7-meson target. A standard analysis of this diagram yields the distribu- 
tion function for heavy quarks (h with mass m) in a heavy meson (H with mass M) in terms 
of a Fourier transform of a heavy quark correlation function along the null plane: 6 

f{x,^) = ^-Jd\e lXx (H(P)\h(0)4h(Xn)\H(P)) fi2 . (2.1) 

Here n M is a null vector defined by = p^ + M 2 n^ /2, the heavy meson momentum and p and 
n are defined by p 2 = n 2 = Q,p-n= 1 yielding p» = p(l, 0, 0, 1), = (l/2p)(l, 0,0, -1). The 
variable x is related to the kinematic variables of Compton scattering by x = Q 2 /2q ■ P. The 
parameter p is arbitrary and fixes the frame. For example, p = Mj 2 corresponds to the heavy 
meson rest frame; p — > oo is the "infinite momentum frame" often used in parton models. 
Here, h(£) is the heavy quark field. Equation (2.1) is written in n ■ A = gauge, avoiding the 
need for a line integral between h and h. Also we have suppressed a factor of n ■ P = 1 in the 
exponent. Contraction with is natural in deep inelastic processes in which some direction 
of hard momentum flow is specified (e.g. = vn^ —xp^ in deep inelastic scattering). It also 
selects operators of pure twist (twist-2 in this case) since traces are removed by contraction 
with n M . The mass scale \i 2 is the renormalization scale required to render the local operators 
in the Taylor expansion of (2.1) finite. Equivalently, /j 2 appears as a cutoff on the transverse 
momentum integral which arises in a momentum space representation of /(x,// 2 ). 
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In order to apply heavy quark symmetry we should take [i 2 = [Xq, where (j,q is a scale 
beneath the lightest of the heavy quark masses, but above the QCD scale so that perturbative 
QCD is still reliable. However, we will see later that to the order of accuracy we can hope to 
achieve, our results apply at any renormalization scale below m 2 . Thus we will eventually set 
Hq = m 2 , but for now we keep it as an independent renormalization scale. The scale m 2 is 
time-like whereas [i 2 is space-like. This distinction only affects terms higher order in a s (/i 2 ) 
than concern us here, so we drop it. 

It is straightforward to convert f(x, (j,q) to a form in which its interpretation as a parton 
momentum distribution is more readily apparent. First, decompose h with the light-cone 
projection operator P± = |7 ± 7 =F with h + = P+h, and insert a complete set of intermediate 
states {|x)}, 

f(x, 1%) = 6 ( x - 1 + P * • n ) l<* \h + (0)\H(P))\ 2 

x 

where we have suppressed reference to the renormalization scale on the right-hand side. The 
momentum P^ = P^ — P£ is the momentum of the struck heavy quark. The function 
f(x) measures the probability of finding a "good" light-cone component of the quark field 
with momentum fraction x (p + = xP + ) in the target. The physical process eH — * eX also 
receives contributions from the heavy antiquark distribution in H, /(x 5/ u 2 ), corresponding 
to Figure 1 with the quark line replaced by an antiquark line. It is difficult to separate / 
from / experimentally. The function f(x,fi 2 ) is defined by the crossing relation f^x.y 2 ) = 
—f(—x,fi 2 ) and has a parton representation like (2.2) with h + replaced by h+. It measures 
the probability of finding a heavy antiquark in the meson H. We argue below that f(x, /a 2 ) 
is small at /s 2 ~ m 2 . Similarly, a potential gluon contribution is also small at [i 2 ~ m 2 . 

-5- 



We now proceed to parameterize (2.1). We first expand 
1 °° 1 r 

f(x,i4) = ^Ey d\e iXx (-i\) k (H(P) |fc(0)n a7 o (n • U) k h(0) \ H(P)) (2.3) 
71 k=o ' J 

where 11^ = iD^ and we have returned to a manifestly gauge invariant form. We match this 
matrix element of the full theory onto matrix elements in the heavy quark effective theory. 
We define the meson four velocity, t> M , by P M = Mw M , where M is the meson mass. Following 
the standard procedure, (although we do not bother to rescale the fields, since everything will 
be dimensionless in the end) we decompose the field h(£) as 

MO = M0 + £*(0 ( 2 - 4 ) 

where 

h v (0 = e-* mv <V + h(0 (2.5) 

h v (0=e— tV-Kt) (2.6) 

where V± = (1± ^)/2. Then 

n M /i(0 = e"^(m^ + n M )(P + /i(0+^-MO) • (2-7) 
We now match the matrix element which appears in (2.3), 

T fc = (H(P) \h(0) ncel a (n ■ U) k h(0)\H(P)) (2.8) 

onto a sum of matrix elements of the heavy quark effective theory. Clearly, these matrix 
elements can be organized as a power series in 1/m. The leading term is the term in which 
the momentum is miy and the matrix element is taken between the fields h v . Mass suppression 
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factors in the remaining matrix elements derive from several sources. First, in the derivative 
in (2.7), the first term is mv^h v (^) = while the second (mass suppressed) term 

is Il^h v (^) w eP^h v {^) where e » 1 — ^ = j^. Second, because the equation of motion 
yields h v (£) = (l/2m)II Q 7 Q /i„(£), the matrix elements involving h v (£) are mass suppressed 
relative to matrix elements of the same operator with h v (£) instead. Finally, matrix elements 
involving the time-ordered product of the operator with mass suppressed operators, Oi/ m , in 
the heavy quark Lagrangian also yield mass suppression. With these considerations in mind, 
we substitute (2.7) into (2.8) and organize the expansion according to the powers of 1/m, 

Tk = Elf \ ) ( n • mv ^ \hv(U)nal a (n ■ U) k ~ £ h v (0)\ H(P)) 



+ U+i) ( n - mv ) t+1 



(H(P) h v (0)n a >y a (n-n) k - e - 1 h v (0) H(P)j 



(2.9) 



+ {H(P) h v (0)n al a (n-U) k - e - 1 h v (0) H(P) 



+ (H(P) T [O 1/m h v {0)n al a {n-Iif- 1 - 1 h v (0)) H(P)) 



+ 



where ... represents terms multiplying ^ , n > 1. In (2.9), (^j^J * s the binomial 
coefficient defined to be zero for £ < or £ > k. The range of the sum on £ is then dictated 
by the values for which {^j^j 7^ 0- The first term has the largest binomial coefficient and 
comes from terms in which all the mass suppression factors come from the small momentum 
piece of II M . In the second term, one power of mass suppression has a different origin. The 
. . . represents terms with even smaller binomial coefficient, which begin at order 1/m 2 . 
We now rewrite the above as 

M k = 2j2( k i )(n-mv) i e k - e A k _ e 

+ ^ k +1 ){n.mvf +1 e k -'B k _, (2.10) 
+ ( ^ 2 Wm^ + V-^ fc _, + ... 



where 



(H(P) h v (0)n al a (n ■ U) k h v (0) H(P)) = 2A k e k (2.11) 



and similarly for B and C, defined in terms of the sum of matrix elements multiplying a 
particular binomial coefficient. * All dependence on the heavy quark mass is contained in the 
factors e p provided e scales inversely to the heavy meson mass. We make the specific choice 
e = A/M = l — m/M which simplifies the algebra. Notice that e scales like the inverse heavy 
meson mass because n ■ P = 1. To get the heavy quark mass expansion in the end requires 
expanding the meson mass M in terms of the quark mass m. Furthermore, because the heavy 
quark occurs in the functions B, C, . . ., one must be careful to expand the mass factor here in 
terms of the meson mass for consistency. With only the functions A, B, and C, the expansion 
above is sufficient. 

Having extracted the factors of e, the parameters A, B and C are independent of the 
heavy quark mass. Here, A = 1 because h rji h counts heavy quarks minus antiquarks in the 
if-meson. 

At this point, we wish to exploit the fact that the same matrix elements appear in all 
the Tfc. We will first add the terms into a general function incorporating this fact. Following 
this, we will derive relations among moments. It is useful to define functions a(y), b(y), etc. 
whose moments give the A k , Bj~, etc. We define 

.4, = / dya(y)y k (2.12) 

oo 

1 

k 



-f 

J — c 

B k+1 = [ dyb(y)y k . (2.13) 
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* Although order e terms can be shown to vanish, this will not be the case for the frag- 
mentation function. We simply leave in all the terms here. 



The limits on the y-integration are chosen to reproduce the required support of the structure 



function (0 < x < 1). Note that a{y) and b(y) must fall rapidly as y — > — oo in order that 
(2.12) and (2.13) converge. The exact nature of this requirement and its physical significance 
will be further discussed below. The upper limit is unity due to our choice of e. So long as e 
scales like 1/M the upper limit is mass independent. We now have 



(2.14) 



// 777 \ ^ / 771 \ ^ 

where we have used n ■ mv = m/M. We can now substitute Tk into (2.3) to obtain 

= £ /~ ^ / 1 oo * e,A ^- a ) fc £ + e? /) fc («(!/) + ^) + ...) 

= — / dy / d\ e ^-m/M-ey) (ffl(y) + c6(y) + _ } 
J — oo J — oo 

= ! S(x- m/M - ey) (a(y) + eb(y) + . . .) 

J — oo 



1 

-a 

e 



x 



777 

M 



x 



+ b 



777 

M 



+ ... 



(2.15) 

(2.16) 
(2.17) 

(2.18) 



Equation (2.18) is the result from which our conclusions about the heavy quark distri- 
bution function follow. It expresses the fact that f(x, Hq) can be expanded in a power series 
in e = A/M. At any order in e, there is a single function. However, the argument of that 
function, 

777 777 
X — X — — 

(2.19) 



y = 



777 ' 



M 



depends in a rather singular way upon e. 



The derivation of (2.18) suggests that it holds only if the functions a(y), b(y), etc. fall 
faster than any power of y as y — > — oo as indicated by (2.12) and (2.13). In fact this 
assumption is stronger than is required. It was convenient for us to expand the distribution 
function in terms of local operators and then resum into f(x). Had we instead matched 
directly onto the bilocal operator product of (2.1), we would have achieved (2.18) directly 
without any requirement that all of the Ak, -Bfc, etc. exist. The steps in such a derivation 
are straightforward, so we summarize them very briefly here. First, perform the matching 
summarized by (2.5) and (2.6) directly in (2.1). The result is 



/ ( Xj /J) = _L J d\ e lXx (H(P) \V+h(0) i V+h(Xn) \ H(P)) e" imAu ' n + . . . , (2.20) 

where we have kept only the term which will eventually turn into a(y) and . . . denotes all other 
terms. Eq.(2.20) remains in n ■ A = gauge (but can be made gauge invariant by restoing 
the line integral) and we suppress the renormalization scale fi 2 . The matrix element in (2.20) 
is dimensionless, remembering the dimension of n^, so converting to the normalization of the 
heavy quark theory it turns into a dimensionless function of the variable A-^j, 

(H(P)\V^(0)4V + h(Xn)\H(P))=f(x^ , (2.21) 

which we represent as a Fourier transform, 

f^\A^ =2 J dae~ iaX ^a{a) . (2.22) 

It is then a simple matter to substitute this Fourier representation into (2.20) and obtain 
(2.18). It is clear from this derivation that we do not yet need to make any assumptions 
about the behavior of a(y) as y — > ±oo in order to obtain (2.18). 



With this second derivation of (2.18) in hand we can examine the behavior of //q) 
near the kinematic limits x = and x = 1. The distribution function /(x,/^) must vanish 
for x > 1 since there are no physical intermediate states in this range. The natural constraint 
o(l) = 6(1) = . . . suffices to guarantee that f(x, /Jq) vanishes at x = 1 independent of e. 

The function f(x, ^q) need not vanish for x < - for negative values of x, it is related 
to /(— x,/j,q). However, the limit x — > corresponds to j/ — > — 1/e in (2.19) and this is a 
physically interesting limit from the perspective of the heavy quark expansion. By comparing 
the arguments of the 5-functions in (2.17) and (2.2), we see that x — > corresponds to 
spectators (light quarks and glue) with momentum P£ ~ m/, since P x • n = e(l — y). Since 
specatators with momenta of order mv^ are not allowed within the framework of the heavy 
quark effective theory, we cannot study the x — > limit in the context of the heavy quark 
expansion. However, we can conclude on the physical basis of the heavy quark expansion that 
the functions a(y), b(y), etc. vanish rapidly as y — > — oo. Any more detailed statement is 
beyond the scope of the heavy quark effective theory. 

In the following section, we will describe a formulation of the heavy quark expansion in 
which the parameters A k , B k , etc. play a central role. The assumption that these coefficients 
exist up to some index k max requires that the functions a(y), b(y), etc. vanish faster than 
y~ k ~ x , y~ k , etc., respectively, as y — > — oo. If a(y), b(y), etc. vanish faster than any power of 
y in this limit, then of course, the A k , etc. exist for all k. If they vanish only like some 
power, then the parameters Af~, B k , etc. exist up to some k max . We cannot determine k max 
from the heavy quark expansion; we take it to be a parameter of our analysis on which we 
comment when appropriate. 

Since all mass dependent factors have been scaled out of A k: B k etc., these should all 
be of order unity. In fact, in free field theory all the A k would equal unity and a(y) would 



be a 5-function at y = 1, corresponding to P£ = 0, in (2.2) (i.e. no spectators). In the real 
world the heavy quark is not free and a(y) is not a 5-function. However, since a(l) = and 
a(y) — > as y — > — oo, it is clear that (2.18) describes a function, /(x,/Xq), which is sharply 
peaked near a; = m/M rs 1 — e and whose width is of order e. In fact, since Ao = 1, /(x, /ig) 
approaches a 5-function as e — > 0, reproducing free field theory in this limit as it should. It is 
important to appreciate the fact that a(y) provides a heavy-quark-mass independent measure 
of the influence of confinement on the heavy quark momentum distribution in a heavy meson 
- a fact which could provide tests of the heavy quark expansion. This will be discussed 
below. 

Equation (2.18) gives f(x, Hq) at the heavy quark mass scale. Experiments will generally 
be done at a value of momentum squared, Q 2 , higher than this. What we have computed so 
far is simply the "boundary data" for QCD evolution from a low mass scale to the scale Q 2 . 
This evolution is non-trivial and complicates the comparison of (2.18) with experiment. We 
postpone the important discussion of comparison with experiment until we have considered 
the fragmentation function. 

II. 2 Fragmentation Function 

The analysis is very similar to that of the previous section. We first define the heavy quark 
fragmentation function as the leading twist contribution of the heavy quark to the process 
e + e~ —* H(P) + X where H(P) is a stable meson with four-momentum P^ containing the 
heavy quark h and X is any additional hadronic contribution to the final state. We first give 
the analog of (2.1) for the spin averaged fragmentation function. 7 

f(x^ 2 ) = ±J2 / d\e^ (0\n al rMXn)\H(P) X ) (H(P) X \h t (0)\0) (2.23) 

x 
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where we have explicitly incorporated the Lorentz indices. As before, we suppress fj 2 on the 
right-hand side and assume n ■ A = gauge. Here, because we cannot perform the sum on x, 
the relevant product of matrix elements does not reduce to a bilocal operator. As with the 
structure function, the process e + e~ — > H(P) + X also receives a contribution from f(x, jj?) 
which measures anti-/i-quark fragmentation into H plus anything. This should be negligible 
for \i 2 ~ m 2 . For now, we retain the variable x = Q 2 /2P ■ Q (1 < x < oo), which in the 
center of mass frame is the ratio of the beam energy to the energy of the hadron. With this 
variable, f(x,Q 2 ) can be interpreted as the probability of finding a hadron of momentum 
fraction 1/x when a quark of unit momentum is created from the vacuum by a current of 
virtuality Q 2 . As before, this can be seen by translating h(Xn) to h(0) and using light cone 
projection operators, 

2 . (2.24) 
* x 

The evaluation of f(x, y?) proceeds similarly to that of f(x, fi 2 ). Here we define 

£ (0 \n al a (n ■ U) k h v (0)\H(P) X ) (H(P) X |^(0)| 0) = 2e k A k (2.25) 

x 

and similarly for the subleading terms. We also define functions a(y) and b(y) by 

/oo poo 
dyy k a(y) , B k+1 = J dyy k b{y) , (2.26) 

etc. 

The range of a(y) has been chosen so that the fragmentation function will only be non- 
zero in the physical region, 1 < x < oo. Following an identical analysis to that for the 
structure function, we derive 



f(x,/4) = 




(2.27) 
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where b is the subleading term defined analogously to the previous section and . . . refers to 
further suppressed functions. With the more conventional definition of the fragmentation 
function variable, z = 1/x, one obtains 



iA 2, 1. 

/(-»^o) = ~ a 

z e 



(I 


m \ 




(I 


m \ 


z 


" M 


+ b 


z 


~ M 




e 






e 


V 


J 




{ 


/ 



+ ... . 



(2.28) 



Here z is the fraction of the beam energy carried by the hadron. 

Equation (2.28) summarizes the predictions of the heavy quark expansion for the frag- 
mentation function at the heavy quark mass scale Hq. As in the case of the distribution 
function, the derivation seems to require the existence of all the parameters Ak, Bk, etc.. 
However, as before, there is an alternate derivation which deals directly with the (2.23) with- 
out ever expanding in local operators. The result is the same, (2.28), without the restriction 
that the integrals defined in (2.26) exist. All of the remarks made in the previous section 
concerning the i ^ or «/ ^ -oo limit apply as well to the fragmentation function in the 
z — > or y — > oo limit. We shall generally assume that the parameters A k , B k , etc. exist for 
k up to some k max . To fit the parameters, (2.28) must be evolved to a scale Q 2 > /j,q. We 
therefore discuss QCD evolution in the next section. 
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III. EVOLUTION, LIMITATIONS OF THE HEAVY QUARK EXPANSION 
AND COMPARISON WITH EXPERIMENT 

IIL1 Moment Analysis 

We have obtained the general form for the fragmentation function at a low renormaliza- 
tion scale. Because the fragmentation function is measured at Q 2 3> m 2 , it is necesary to 
discuss the QCD evolution between the two scales. This is most conveniently accomplished by 
evolving the measured moments of the fragmentation function to the scale Hq at which we have 
parameterized the function. At this scale, one extracts the mass independent parameters. 

In this subsection, we will discuss the application of this procedure to the first few 
moments of the fragmentation function. In later subsections, we explain why we restrict our 
analysis to only these lowest moments. We parameterize these moments in terms of mass 
independent parameters (related to those which entered the functions a(y) and b(y)) and 
show furthermore how to extract the parameters charactizing the moments at leading and 
sub leading order in A/ra. 

We define the moments of the fragmentation function at a scale Q 2 above the heavy 
quark scale in the usual way. 

fk(Q 2 )= f dzz k f(-,Q 2 ) . (3.1) 
Jo z 

Notice we study the moments in z for positive k. These are the conventional moments: they 
scale simply under QCD evolution and are convergent at all scales Q 2 . Having defined the 
moments at the scale Q 2 , one must renormalize them from Q 2 to the scale /z 2 ,, where we have 
parameterized the matrix elements. This proceeds in three steps. First, one must evolve the 
function from the high momentum scale Q 2 at which the experiment takes place to the heavy 
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quark mass scale m 2 . Then one must match the evolved moments onto the heavy quark 
theory. Finally, one must evolve in the heavy quark theory from the scale m 2 to the scale \i\ 
at which the matrix elements were parameterized. To apply heavy quark flavor symmetry, 
this scale should be the same for the b and c quarks. 

The first stage of evolution, between Q 2 and m 2 , is standard. The fragmentation function 
measured in e + e~ — > H + X at the large scale Q 2 is a linear combination of quark, antiquark 
and gluon fragmentation functions. The gluon and antiquark fragmentation functions are not 
zero at the scale Q 2 . However, they evolved from the scale m 2 at which both the antiquark and 
gluon fragmentation functions were negligible. In particular, because the quark is heavy, the 
gluon moments at the scale m 2 are suppressed by (A/m) 2 relative to the quark contribution. 
With these constraints on the antiquark and gluon fragmentation functions it is possible to 
de-evolve the data from Q 2 using the standard anomalous dimension matrix and obtain the 
heavy quark fragmentation function at the scale m 2 . Above the heavy quark mass scale, 
the fact that the quark is heavy is irrelevant to the evolution equations, and enters only in 
the boundary condition which tells us the gluon moments are small at the scale m 2 . This 
information about the boundary data, together with the experimental data at the scale Q 2 , 
is sufficient to solve completely the evolution equations. 

The second stage in the procedure is to match from the full theory onto the heavy quark 
theory. From the definition of the fragmentation function given in (2.23), it is clear that 
the moments match onto the moments in the heavy quark theory in a straightforward way. 
That is, we already know how to match the matrix elements appearing in the function itself 
onto matrix elements of the heavy quark effective theory characterized by the parameters 
which we have defined. The moments are defined as integrals of f(x,Ho) times 



powers of x. This definition holds in both the full and effective theories. If one were to work 
to subleading order, there would be additional contributions from evaluating the one loop 
matrix elements in the full theory. 

The final stage is to de-evolve the moments in the heavy quark effective theory from the 
heavy quark mass scale, m 2 , to the uniform scale of the effective theory, /1q. Fortunately, to 
leading and sub-leading order in e, this is trivial. The point is that the matrix elements at 
leading and subleading order in the heavy quark mass expansion do not evolve. This is most 
easily seen in v ■ A = gauge where is the heavy meson velocity. Because the parameters 
Ak and Bk are defined in terms of matrix elements of local, gauge-invariant operators in 
the effective theory which involve only a single velocity v^, the gluons do not couple in 
leading order in the heavy quark mass expansion. Of course, at the level of 1/m 2 , there are 
QCD radiative corrections in the heavy quark effective theory. Therefore, in order to extract 
parameters at order e 2 further QCD evolution would be necessary. 

To summarize, the moments of the heavy quark fragmentation function should be de- 
termined experimentally from the data at the scale Q 2 . They should then be de-evolved via 
QCD evolution down to the scale m 2 , at which point they should be matched to the mo- 
ments of the function f(x,m 2 ) parameterized in the heavy quark effective theory through 
the parameters Ak and B\, or equivalently, the functions a(y) and b(y) which summarize the 
relations among the moments. The QCD radiative corrections in the full theory are given 
in Ref. [5] to subleading order. For the leading parameters no further QCD evolution in the 
heavy quark effective theory is required. 

For simplicity, we therefore compare the evolved, experimentally determined moments 
at \j? = m 2 to the moments derived from the function (2.23). It is more convenient to work 



with the variable x to which we now revert. Then 



/oo 
dxx~ k ~ 2 f(x,m 2 ) 



m \ / m 



oo / 2 / X — — \ /X 



efez-*- 2 -a ^ +6 ^ I + ... 



(3.4) 



e 



A simple change of variables yields 

dy [ey + — J (a(y) + e6(y) + . . .) 
where we have shown explicitly only those functions which begin at order e° or order e. Here 



(3.5) 



oo 



/ dy{y-lfa{y) 

J \oo (3-6) 
B m = y dyfo - ifbiy) . 

Notice that only a finite number of Ag, Be, etc. occur in the expansion of any given moment 
to a fixed order in e. 

Consider now the explicit expressions for the first few moments. 

f (m 2 ) = A + e (-2A 1 + B t ) + e 2 (3l 2 - 2B 2 + C 2 ) + . . . 

f x(m 2 ) = i + e (-3ii + Si) + e 2 (6i 2 - 3B 2 + C 2 \ + . . . 

(3.7) 

f 2 (m 2 ) = A + e (-41 1 + Si) + e 2 (l0i 2 - 45 2 + C 2 ) + . . . 

f 3 (m 2 ) = i + e (-5ii + Si) + e 2 (l5i 2 - 5S 2 + C 2 ) + . . . 

If we work to order e and require the remainders in (3.7) to be bounded, then it is necessary 
to assume that a(y) vanishes faster than y~ 3 for large y. Likewise b(y) must vanish faster 
than y~ 2 . In order to work to order e 2 , a(y) and b(y) must vanish one power faster for large 
y. Assuming these asymptotic forms, and that the heavy quark expansion is valid, we can 



extract the leading parameters of the functions a(y) and b(y). In particular, from knowledge 
of either the b or c quark fragmentation function alone, one can extract the parameters 
A Q + eBi + e 2 C 2 + . . . ~ A , eA\ — e 2 B 2 + ...?» eAi, and e 2 A 2 , working only to order e 2 . 
With knowledge of both the c and b quark fragmentation functions, one can independently 
extract Aq, eBi, eAi, e 2 B 2 , and e 2 A 2 . Of course if we work only at order e, there is an equal 
spacing rule which serves as a consistency check on the moments: 

f 1 -f = f 2 -f 1 = f 3 -f 2 . (3.8) 

We can now also see the main problem with extracting information from higher moments. 
In order to extract a parameter at any specified order in the heavy quark mass expansion, 
we need to assume that higher order terms can be neglected. However, the large binomial 
coefficients multiplying higher order terms in e make this impossible for the higher moments. 
That is, for the fcth moment, with k large, the expansion is not in e, but in ke, which 
compromises the application of the expansion beyond the lowest moments. In the next sub- 
section, we discuss the further advantages to be had by restricting attention to only the 
lowest moments. After that (in III. 4), we return to the question of whether information can 
be extracted from the higher moments, or equivalently from the function f(x,m 2 ) itself. 
III.2 Subleading QCD Evolution 

We have discussed leading logarithmic QCD evolution. However, it is well known that 
the usual QCD evolution program has difficulties with heavy quark fragmentation functions. 
In this section, we explain why these problems are less severe for our program than for a fully 
perturbative study of heavy quark fragmentation functions and may be ignored entirely if we 
restrict our analysis to the lowest moments. 
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In principle, one would expect that QCD evolution in the full theory between Q 2 and 
m 2 is perturbative, as both scales are much higher than the QCD scale. However, near 
x = 1, the QCD coupling a s is multiplied by logarithms of the form log(l — x), etc. In this 
region, effects which are formally subleading are nevertheless large. In fact, it has been shown 
that the summation of these Sudakov logarithms can be summarized by evaluating a s at the 
scale (1 — x)Q 2 . 8 These effects are important when (1 — x)Q 2 « Q 2 and furthermore the 
analysis suggests that other, truly non-perturbative effects cannot be neglected for x such 
that (1 - x)Q 2 ~ A 2 . 

Now consider the form of the function renormalized at m 2 given by (2.28). It vanishes 
at x = 1 and is peaked at x — 1 ~ e. Sudakov effects are worst at the lowest Q 2 , in this 
case Q 2 = m 2 . At Q 2 = m 2 , non-perturbative effects beyond the resummation of Sudakov 
logarithms are important for 1 — x < e 2 . The distinction between e and e 2 is crucial. Since / 
is likely to be a rapidly falling function in the region close to x = 1, and since the scale for 
variation in / is e, the function should be very small in the region where Sudakov logarithms 
are important. In particular, the first few moments are very likely dominated by x in the 
regime where even Sudakov logarithms are not large, so a leading logarithmic analysis should 
be valid. However, we do not know how fast / falls as i — > 1. This non-perturbative 
information determines which moments are trustworthy. Since constraints on the heavy quark 
expansion most likely restricts us more severely to low moments than do considerations of 
Sudakov effects, we expect that Sudakov logs can be neglected altogether. 

This is in contrast to previous works 4 ' 5 which attempted a fully perturbative analysis 
of the fragmentation function over the whole range of x. They have been led to a form of 
the fragmentation function which is singular at x = 1 and have necessarily had to deal with 
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evolution in a regime in which non-perturbative effects are important. From this perspective, 
our approach has been to summarize the nonperturbative information in a few parameters. 
With sufficiently rough averaging (by looking only at low moments) we can then do a simpler 
QCD analysis. 

We see now the essential advantage gained by restricting attention only to the lowest 
moments, independent of the breakdown of the e expansion. Our objective is to cleanly 
distinguish perturbative and non-perturbative QCD. For x very close to 1, this distinction is 
not possible. But this region of x should not contribute significantly to the low moments for 
the reasons outlined above. 

Having addressed the issue of Sudakov logarithms, the study of higher order QCD effects 
should be standard. For example, there is an ambiguity associated with the precise scale at 
which we match between the full and heavy quark effective theories which can only be resolved 
through a higher order QCD calculation. Moreover, as a matter of practice, the data on the c 
quark is obtained at center of mass energy of the T(4S) Furthermore, it might be that in this 
case the logarithm is not sufficiently large for a leading log analysis to be adequate. Certainly, 
a sub leading calculation would be required to extract a parameter at the precision of e 2 . 

III. 4 Direct Analysis of f(x 7 m 2 ) 

Having analysed the limitations on the moment analysis, we return briefly to the pos- 
sibility of extracting information from looking directly at the function f(x,m 2 ). As long as 
we restrict (1 — x) 5>> e 2 Sudakov effects are not likely to be important. Looking at (3.7) we 
see that the moment analysis breaks down when ke > 1. However, the terms which become 
large in this limit can be resummed yielding the original functional form, (2.28), in which 



successive terms are suppressed by powers of e when the argument, y = — is held 

fixed. So it is appropriate to analyze f(x, m 2 ) directly. 

First, suppose the function f(x, m 2 ) can be extracted from experiment and QCD de- 
evolution. The predictive power of (2.28) lies in the orderly expansion in e implied for f(x, m 2 ) 
when x = 1 + e(y — 1), 

/ (1 + e(y - 1), m 2 ) = h(y) + b(y) + ... (3.9) 

where a(y), b(y), etc. are all independent of e. In principle, / can be measured for both c- 
and 6-quarks, a(y), b(y) and e can all be extracted and the heavy quark expansion can be 
tested. To leading order in e, (2.28) requires that there is a choice of A (which in principle 
could have been extracted elsewhere 10 ) such that 

ej c (l + e c (y-l),m 2 c ) - e b f b (l + e b {y -l),mg) =0 (3.10) 

Having fixed e c and e b , a(y) can be extracted from 

a(y) = (f h (1 + c 6 (l - y), m 2 ) - f c (l + c c (l - y), m 2 )) (3.11) 

A cautionary note is appropriate: After de-evolution f(x, m 2 ) would likely be available only 
with some finite resolution (smearing) considerably greater than the original experimental 
resolution at Q 2 . If the resolution were as large as the intrinsic width of f(x, m 2 ), which is of 
order e, then it is easy to see that this method loses its usefulness and the moment method 
outlined previously gives as much information as can be extracted from the data. 

III. 3 Discussion 

There are several points we have neglected in the previous analysis which we now clar- 
ify. First, we have made assumptions about the behavior of the heavy quark and antiquark 
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distributions at the heavy quark mass scale. In particular, we assumed f(x,m 2 ) vanishes 
rapidly as x — > oo and we have ignored f(x,m 2 ). These assumptions are only sensible if 
they are stable under perturbative QCD evolution in the heavy quark theory. They would 
clearly not be consistent above the heavy quark scale, at which point conventional evolution 
generates a tail in / which reaches to z = 0, a gluon distribution at order a s ln(<5 2 /m 2 ), and 
a heavy antiquark distribution / at order a 2 ln 2 (<5 2 /m 2 ). Thus it would not be reasonable 
to implement such assumptions in the domain of perturbative evolution. However, we are 
placing these restrictions on the boundary data at and below the heavy quark mass scale, m 2 . 
As we have already discussed, the structure function does not evolve at scales below m 2 at 
leading or subleading order in the heavy quark mass expansion. At this level of accuracy, our 
assumptions are therefore stable with respect to with perturbative QCD evolution. 

We have also neglected gluon and antiquark contributions to the moments at the scale 
m 2 . Because a gauge invariant gluon operator matrix element is suppressed by at least e 2 and 
receives no enhancement by a binomial coefficient relative to the leading and subleading terms, 
it can only contribute to the function c(y) and other suppressed functions. At the level in the 
heavy quark expansion at which we are working, we can neglect the gluon moments at the 
scale m 2 . However, as we have discussed, gluon moments are generated through perturbative 
QCD evolution so they do not vanish at the scale Q 2 . Similarly, antiquark moments are only 
generated through QCD evolution. 

The next point to clarify is that throughout our analysis, we have assumed fragmentation 
into a specific final state, namely the pseudoscalar meson. However, experiments at very high 
center of mass energy will not distinguish the pseudoscalar and vector B mesons. Further- 
more, the final state could be a B s meson or a baryon containing a b quark. However, the 
fragmentation function into D* is what is best measured at CESR. 
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If we are to apply our methods to compare the fragmentation function at LEP and 
CESR, we would not expect to obtain results which are accurate to better than order 20%. 
This is because heavy quark symmetry only guarantees the identity of the matrix elements 
between the pseudoscalar or the vector states at zeroeth order in 1/m. This means that the 
parameter eB\ is not identical for the pseudoscalar and vector mesons. However, counting 
spins shows that the vector states are produced three times as often as the pseudoscalar. So 
even if the time ordered product of the spin splitting operator and the zeroth order matrix 
element accounted for 50% of the contribution to eB\, the relations between the B system 
and the D system (namely the vector states) should be accurate at the 15 to 20 % level, 
which is the inaccuracy represented by higher order terms in the heavy quark expansion. 
(The mass splitting between the pseudoscalar and meson masses is higher order, and can be 
safely neglected.) 

B s production should also not be a problem at this order of accuracy. The B s is only 
produced about 12 % of the time, and chiral symmetry should guarantee that the matrix 
elements are equal to those of the Bj at the 20% level. 

Finally, baryons will probably be produced about 10% of the time. This can be a problem 
since the parameters of the baryons are independent of those of the mesons. The danger is 
that an inaccurate determination of Aq would lead to a large error in extracting eA\ and eB\. 
However, there is only a single new parameter at the level of accuracy of our calculation, 
namely the fraction of baryons, p (since pe is already only a few percent.). If a guess of p at 
about 10% is accurate at the few percent level, one should still be able to extract the same 
parameters as before. If p is left a completely free parameter, eB\ can still be extracted. And 
of course if the baryons can be excluded from the sample, everything would proceed as in the 
text. In fact, the same analysis would apply to baryons as to mesons. 



Furthermore, if it is possible to study D mesons without restricting oneself to _D*'s, but 
instead sums inclusively over all states, the relative fraction of mesons and baryons should 
be independent of heavy quark mass at leading order in the heavy quark mass expansion, so 
that the parameter p need not be extracted. The analysis we gave would apply, but with each 
parameter in the moment expansion interpreted as a sum of a meson and baryon contribution. 

It should be kept in mind that the consistency tests among the moments of any particular 
flavor quark should hold even with different final state particles, since they only required the 
existence of a heavy quark expansion, which is true for each final states individually. It is 
only when comparing between different heavy quark systems with different fractions of the 
various final states that complications arise, unless one can identify the contributions from 
the particular final states. 

Notice also that we have expanded in e, which scales inversely to the meson mass so the 
expansion is not strictly in heavy quark mass. It is straightforward to convert our results to 
an expansion in heavy quark mass. 

Of course, there are potential experimental difficulties which we expect could make an 
analysis beyond the level of e impracticable. In particular, finite experimental resolution 
could be a problem. The spread of the measured fragmentation function of the b is increased 
significantly by the perturbative QCD contribution. The measurement has to be sufficiently 
accurate to determine the nonperturbative parameters once the perturbative contribution has 
been removed. Furthermore, what is measured is the hard lepton in the final state. With no 
further theoretical input one needs an accurate relationship between the momentum of this 
lepton and the momentum of the meson itself in order to accurately measure the fragmentation 
function. This requires good measurements of the decay lepton spectrum from B decay. 



Alternatively, one would require a large sample of fully reconstructed events. However with the 
heavy quark expansion one could determine the lepton energy distribution, thereby simplifying 
the analysis. 

Finally, we comment on the commonly used Peterson fragmentation function. It is clear 
from any analysis of heavy quark fragmentation functions that one expects a function which 
is weighted towards 1 but with some spread, and which goes to zero at z = and z = 1. 
Any function of this sort will more or less fit the data. So the Peterson function is as good 
as any other. However, the Peterson function assigns a width of order e 2 , where we have 
seen the true width of the fragmentation function is order e. Therefore, one must be careful 
in assigning a physical interpretation to the parameters of the fragmentation function, or in 
trying to relate the fragmentation functions for different flavors of quark. 
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IV. CALCULABLE HIGHER-TWIST CORRECTIONS 

We now turn to the subject of higher-twist corrections to e + e~ — > H(P) + X. In general, 
the study of higher- twist effects is theoretically complicated and experimentally difficult. In 
this case, however, the heavy quark mass sets the scale for higher-twist effects which are poten- 
tially the most important because they are of order m 2 /Q 2 which in certain experiments can 
be large (e.g. m b = 5 GeV, Q = 15 GeV gives m 2 /Q 2 > 10%). Fortunately they can be related 
to the leading twist fragmentation function. Because the subject is technically complicated, 
we first summarize our results for the reader who does not wish to see the derivations. We 
next give a simple, non-technical derivation of the primary results. Finally we summarize the 
results of a complete, operator-product expansion through twist-four (0(l/Q 2 )) distribution 
functions which will be published elsewhere. 10 

First, our conclusions: At leading order in a s , we show that there are no power corrections 
to e + e~ — * H(P) + X which scale by the heavy quark mass alone, that is, no corrections of 
the form (m 2 /Q 2 ) k for k > 1 to either the transverse or longitudinal fragmentation functions 
defined by 

*V=^E / d^e^(0\J,(O\H(P)X)(H(P)X\M0)\0) (4.1) 
( QnQv\ - / 2 \ 

+ [g U v ~ + ^2 — ] w T{q,v), (4.2) 

where v = P • q. [The "fragmentation function" discussed in Section II is the scaling limit of 
wt(q 2 , v). The function wl{<1 2 , v) vanishes at leading twist on account of a Callan-Gross-like 
relation.] At higher order in a s (Q 2 ), there are radiative corrections which yield twist-four 
corrections of the form 

m 2 fa s (Q 2 ) ln Q 2 



Q 2 \ n m 2 



(as well as nonvanishing wl(Q 2 )). Because of heavy quark constraints on matrix elements at 
the heavy quark scale it is possible to calculate the magnitude of these corrections exactly 
in terms of the leading twist heavy quark fragmentation function described in Section II. 
There are also corrections of order mA/Q 2 which can also be related to the leading twist 
fragmentation function and are calculable if A is known. 

We conclude that the process e + e~ — > H(P) + X will probably not have large higher 
twist corrections, even at relatively low center of mass energy. At heavy meson threshold, 
Q 2 = AM 2 , ordinary higher- twist corrections are suppressed by A 2 /Q 2 rs A 2 /AM 2 « 1 so 
these are cleary small. Corrections of order mA/AM 2 are only of order 10 % for the c quark 
and smaller still for the b quark. However, this could confuse the extraction of an order e 2 
parameter. Corrections of order m 2 /AM 2 (a(4M 2 )/-7r) In (4M 2 /m 2 ) are less than 5% for the 
c-quark and much smaller for the 6-quark. These should be well under control. 

Our results follow from simple considerations of the structure of heavy quark mass cor- 
rections to e + e~ — > H(P)X. The analysis is conceptually simpler for the deep inelastic 
scattering process eH(p) — > eX where one can use the language of the operator product ex- 
pansion. We will present our argument for deep inelastic scattering but it should hold equally 
well for e+e~ -> H(P)X. 

We first show that there are no corrections of the form (m 2 /Q 2 ) k . We work to lowest 
order in a s . Since we are interested in effects of order (m 2 /Q 2 ) fc , we may also set A = 0. 
Thus, we may ignore QCD entirely and calculate eH(P) — > eX from a free, on-shell heavy 
quark h as shown in Fig. 1. This calculation, though trivial, gives us the information we 
require. This is true only because the relevant matrix elements are equal to their free field 
values up to non-perturbative corrections suppressed by e. So at order e° and ck°, a free field 
theory analysis suffices. 
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Direct calculation of Fig. 1 gives 

Wi»< = < P )^ d + V +m) S [(q + Pf - m 2 ] lv u{P) (4.5) 

spins 

= - A Tr {If +m) ( 7 „ i +2P M ) 5 (q 2 + 2q ■ P) lv 

If P^q, + Puq^ P„Pvq 2 \ s , ^ 
= g + ^2— J S(x - 1) (4.4) 

which is manifestly independent of the quark mass. Notice that because we have done this 
calculation in free field theory, the structure function itself is trivial: fi(x) = \5{x — 1) in the 
limit a s , A — > 0. This is the structure function of a free quark, the inverse Mellin transform 
of the matrix elements of the only twist-two operators which persist in this extreme limit: 

^i-Mn = «7 Mi n^ 2 . . . II Mri /i - traces (4.5), 

all of which have identical coefficients (A k = 2) in the operator product expansion. The 
matrix elements of the O n will, in reality, deviate from the free field limit by corrections of 
order A/m — already described in Section II — but these modifications will not alter the 
relations among coefficient functions which result in the cancellation of all O (m 2 /Q 2 ) terms 
from (4.4). 

To understand this result more deeply and to see why the corrections of order 
m_ ) } n Q_ j anc | Am/Q 2 are calculable, it is necessary to consider the process from 

the viewpoint of the operator product expansion. There are many sources of Q 2 suppressed 
corrections. The situation is made more tractable by working at tree level, postponing the 
study of the effects of evolution and operator mixing which will of course generate corrections 
of order a s to this result. At tree level and through order 1/Q 2 there has been a complete 
analysis of electroproduction with massless quarks. 11 The twist-four operators which appear 
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in that analysis include four-quark operators, generically (hh) 2 , two quark gluon, hFh, and 
two quark two gluon, hFFh, operators. It is easy to see that, renormalized at a scale where 
the heavy quark effective theory applies, these operators contribute only at order A 2 /Q 2 , The 
dangerous terms are those operators directly proportional to m or m 2 which were omitted 
from the calculation of Ref. [11]. Four categories of such contributions occur: 

1. "Target mass corrections" — These well-known kinematic corrections were first under- 
stood completely by Nachtmann. 12 They arise when twist-two operators — traceless and 
symmetric in all indices — are projected out of the dominant operators which appear 
naturally in the expansion of two currents. Target mass corrections are proportional to 
(M 2 /Q 2 ) fc where M is the heavy meson mass, and their x-dependence is determined by 
the leading twist-two distribution or fragmentation function. 

2. Modifications of the struck quark propagator — The replacement 1/p 2 — > l/(p 2 — m 2 ) 
generates corrections of order (m 2 /Q 2 ) k to the coefficient functions of the leading twist- 
two operators. 

3. New twist-three operators with coefficient functions linear in m — The quark mass term 
in the Dirac propagator, {j> +m)/(p 2 — m 2 ), generates a tower of twist-three, chiral odd 
operators of the form ipD^ . . . D Mn ip with coefficient functions linear in m. The matrix 
elements of these operators are proportional to M and the resulting corrections to the 
distribution and fragmentation functions are of order mM/Q 2 . 

4. When twist-four operators are resolved into a canonical form in which gluon fields are 
manifest, it is necessary to use the QCD equations of motion. The use of the massive 
Dirac equation {If —m)h = generates new twist- four operators proportional to m and 
m . 
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A careful study of the operator product expansion at tree level through twist-four shows these 
are the only corrections through order m 2 . 

The result of the operator product expansion calculation for the distribution function is 
most conveniently expressed as moments of the transverse and longitudinal structure func- 
tions, Ft and Fl, defined by 



F T (x,Q 2 ) = —F 2 (x,Q 2 ) 

F L (x, Q 2 ) = ^ (l + ) F 2 (x, Q 2 ) - F 1 (x 1 Q 2 ) . 



(4.6) 



Then our results are 



Ml 



(m 2 ) = f 
Jo 



dxx k 1 F T (x,m z ) 



T fc (m 2 ) + 



1 fc-2 
Q 2- k~ 



2^ , 2n (k-l)(k+6)« 



(k-3)M 2 T k (m 2 ) - 



k 



m r fc _ 2 (m ) 



f — ; 2mMA fc _ 2 (m 2 ) 



Mi 



(m 2 ) = f 
Jo 



k 



dxx^F^x.Q 2 ) 



(4.7) 



^ ^ |-M 2 r fc (m 2 ) - ^m 2 r fc _ 2 (m 2 ) + 



fc-1 
k 



2mMA k - 2 (m 2 ) 



(4.8) 



where fc = 2, 4, 6, . . . and T k (m 2 ) and Ajt(m 2 ) are the matrix elements of twist-two and twist- 
three heavy quark operators, respectively (at the renormalization scale m 2 ): 



T k (m 2 ) = (H(P) hi (II- nf~ 1 h H(P) 



MA k (m 2 ) = (H{P) h(U-n) k h H{P) 



(4.9) 
(4.10) 



The first term in (4.7) is the standard twist-two result. The other terms in (4.7) and (4.8) 
each consist of three parts. The first, proportional to M 2 I\, are target mass corrections. The 



second, proportional to m 2 I\_2, contain mass corrections to the struck quark propagator 
and operators generated by use of the Dirac equation. The third, proportional to mMkk-2, 
originate in chiral odd, twist-three operators and also from use of the Dirac equation. The 
fact that these corrections can be grouped into transverse and longitudinal contributions is 
a consequence of the gauge invariance of virtual photon if-meson Compton scattering and 
is a highly non-trivial check on the algebra. So is the fact that (4-7) and (4-8) reduce to 
the free quark result, (4-4)? when the {T^} and {A^} are set equal to their free quark values: 
r fc (free) = A fc (free) = 2. 

In general, the {Tfc} and can be related by the QCD equations of motion: 

k-2 

m 2 h i (n ■ Tlf^h = mh(n ■ U) k h + m^%. II) V Q ^F 7/3 n Q n 7 (n • U) k - £ ~ 2 h (4.11) 

£=0 

for k > 1, and 

m 2 h rji h = mhn ■ Tlh 

for k = 1. In general (e.g. for light quarks) the third operator in (4.11) obstructs any attempt 
to relate Tk to Afc. However, at the heavy quark scale it is easy to see that the matrix element 
of the third operator is of order A 2 /m. Define 



C {k) EE ( H(P) 



k-2 



J2 Hn ■ U) i a a f3 F^n a n 1 (n ■ U) k ~ £ - 2 h 



H(P)) . (4.12) 



e=o 

At the heavy quark scale, numerator factors of m arise only from the substitution IP = 
nv^ + plus one additional factor from the norm of the Dirac fields. Each factor of n M 
contributes a factor 1/M. The result is oc m fc_1 /M fc , but on dimensional grounds 
scales like mass, so we conclude 

C'"^' 1 ^ (4.13) 
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to leading order in the heavy quark mass expansion. With this result in hand, we take the 
matrix element of (4.11) in a heavy meson state and find 

m 2 T k (m 2 ) = mMA k (m 2 ) + C(A 2 ) (4.14) 

where we have restored the renormalization scale (m 2 ) as a label on the operator matrix 
elements to remind ourselves that (4.14) holds at the heavy quark scale. Substituting (4.14) 
into (4.7) and (4.8) we obtain the dominant twist-four corrections to the longitudinal and 
transverse moments at the heavy quark scale, m 2 : 



M T k (m 2 ) 
M L k {m 2 ) 



^ {k 2 f 3) [M 2 Y k {m 2 )-m 2 Y k _ 2 {m 2 )]+0{K 2 /Q 2 ) (4.15) 



twist— 4 

4 k- 1 



[M 2 T k (m 2 ) - m 2 r fc _ 2 (m 2 )] + 0(A 2 /Q 2 ) (4.16) 



twist -4 Q k 

Tk(m 2 ) can be expressed in terms of Aq, A\ and B\ (to order e) by means of (3.7). m and 
M are related by m = (1 — e)M. Altogether 

M 2 T k {m 2 ) - m 2 r fc _ 2 (m 2 ) = 2eA x + 0(e 2 ) (4.17) 



so we obtain 



, M 2 „ (jfe-2)(ife-3) , 
Mi(m') =^-^A 1 y - £ i (4.18) 

twist -4 V K 



t L, 2^ M 2 , k—1 



twist -4 V ft ' 



Equations (4.18) and (4.19) hold at the heavy quark scale, Q 2 = m 2 , and at higher scales, 
Q 2 > m 2 , modulo QCD radiative corrections which are discussed below. Like other results 
which make use of the mass expansion, (3.7), they are valid only for low moments, for which 
k <K 1/e because we ignored k 2 e 2 compared to ke in the derivation. 
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Equations (4.18) and (4.19) summarize our assertion that the leading higher-twist correc- 
tions to the scaling of distribution functions is of order AM/Q 2 , and is calculable in terms of 
the parameter eAi, which can, in principle, be extracted from the leading twist heavy quark 
distribution function. From the operator-product-expansion point-of-view, the cancellation 
of terms of order m 2 /Q 2 among their several sources seems fortuitous. The simple graphical 
argument given at the beginning of this section gives a firm physical basis for the cancellation. 
It also assures us that the same result holds for fragmentation functions for which operator 
product expansion methods are less straightforward. 

Finally, it is necessary to comment on the effect of QCD radiative corrections on these 
results. We focus on the moments, Ml(Q 2 ) and M%(Q 2 ), defined by (4.7) and (4.8). The 
operators parameterized by Tk(Q 2 ) and Ak(Q 2 ) have different anomalous dimensions so the 
cancellation of m 2 /Q 2 terms, which was complete at m 2 , is incomplete at Q 2 . The terms 
which survive are clearly O Qs ^ In ^ ) • To establish the way to compute these terms 
it is necessary first to examine the Q 2 values of interest. For Q 2 » m 2 one could use 
standard, leading order renormalization group evolution. However, for Q 2 y> m 2 twist-four 
effects are unlikely to be measurable in the first place. We are more likely to be interested 
in Q 2 >m 2 where twist-four effects might be measured. At such scales there is no reason to 
ignore O (a s (Q 2 )/7i) corrections to coefficient functions. Fortunately, there is a direct way to 

2 q2 

compute all relevant corrections of order ^ ^ In since A may be taken to zero in these 
calculations, they may be extracted from the Feynman diagram calculation of the appropriate 
process involving free, on-shell heavy quarks — the same method we used to show corrections 
of order m 2 /Q 2 vanish at the beginning of this section. The relevant diagrams are shown in 
Fig. 2. The processes shown in Figs. 2b - 2d, which involve mixing of gluonic and quark-gluon 



operators in the language of the operator product expansion, are easily seen to be suppressed 
by A 2 /Q 2 . Thus it is a matter of straightforward calculation to compute the set of all im- 
portant 1 / Q 2 corrections to heavy quark fragmentation or distribution functions: corrections 
of order AM/Q 2 have been calculated in this paper; corrections of order ^p- ^ In can be 
obtained from calculation of Feynman diagrams involving on-shell heavy quarks. 

V. CONCLUSION 

There are only very few places where one can hope to extract mass suppressed contri- 
butions to matrix elements. The existence of these additional parameters should be very 
important from the point of view of testing heavy quark symmetry at subleading order. Fur- 
thermore, it will be useful to see what the size of mass suppressed corrections really turns out 
to be in order to test the validity of the heavy quark expansion. 

Clearly the ideas of this paper need to be tested by application to existent and future 
data. Future work would involve incorporating subleading evolution and calculation of the 
perturbative higher twist effects. 
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FIGURE CAPTIONS 

Figure 1: Virtual Compton scattering of a current coupled only to /i-quarks from a 
heavy meson, H, with momentum P. The imaginary part is proportional to the distribution 
function, f(x,Q 2 ), at large-Q 2 . Only the dominant process at large-Q 2 is shown, exclusive of 
QCD radiative corrections. 

Figure 2: The cross section for an /i-quark current to produce an iif-meson plus anything 
in the deep-inelastic limit. Only the dominant contribution at large-Q 2 is shown, exclusive of 
QCD radiative corrections. The cross section is proportional to f(x, Q 2 ). 
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